Lecture 21
How long is the flight between the Western most and the Eastern most points in the US?
R has a rich package ecosystem for read/writing, manipulating, and analyzing geospatial data.
sp - core classes for handling spatial data, additional utility functions - Deprecated
rgdal - R interface to gdal (Geospatial Data Abstraction Library) for reading and writing spatial data - Deprecated
rgeos - R interface to geos (Geometry Engine Open Source) library for querying and manipulating spatial data. Reading and writing WKT. - Deprecated
sf - Combines the functionality of sp, rgdal, and rgeos into a single package based on tidy simple features.
raster - classes and tools for handling spatial raster data.
stars - Reading, manipulating, writing and plotting spatiotemporal arrays (rasters)
See more - Spatial task view
sfThis is the hardest part of using the sf package, difficulty comes from is dependence on several external libraries (geos, gdal, and proj).
Windows - installing from source works when Rtools is installed (system requirements are downloaded from rwinlib)
MacOS - install dependencies via homebrew: gdal, geos, proj, udunits.
Linux - Install development pacakages for GDAL (>= 2.0.0), GEOS (>= 3.3.0), Proj4 (>= 4.8.0), udunits2 from your package manager of choice.
More specific details are included in the repo README on github.
st_read / st_write - Shapefile, GeoJSON, KML, …
read_sf / write_sf - Same, suports tibbles …
st_as_sfc / st_as_wkt - WKT
st_as_sfc / st_as_binary - WKB
st_as_sfc / as(x, "Spatial") - sp
# A tibble: 4 × 3
path type size
<fs::path> <fct> <fs::b>
1 …a/gis/nc_counties/nc_counties.dbf file 41K
2 …a/gis/nc_counties/nc_counties.prj file 165
3 …a/gis/nc_counties/nc_counties.shp file 1.41M
4 …a/gis/nc_counties/nc_counties.shx file 900
Reading layer `nc_counties' from data source
`/Users/rundel/Desktop/Sta523-Fa22/website/static/slides/data/gis/nc_counties'
using driver `ESRI Shapefile'
Simple feature collection with 100 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32186 ymin: 33.84175 xmax: -75.46003 ymax: 36.58815
Geodetic CRS: NAD83
Simple feature collection with 100 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32186 ymin: 33.84175 xmax: -75.46003 ymax: 36.58815
Geodetic CRS: NAD83
First 10 features:
AREA PERIMETER COUNTYP010 STATE
1 0.11175964 1.610396 1994 NC
2 0.06159483 1.354829 1996 NC
3 0.14023009 1.769388 1998 NC
4 0.08912401 1.425249 1999 NC
5 0.06865730 4.428217 2000 NC
6 0.11859434 1.404309 2001 NC
7 0.06259671 2.106357 2002 NC
8 0.11542955 1.462524 2003 NC
9 0.14328609 2.397293 2004 NC
10 0.09245561 1.810778 2005 NC
COUNTY FIPS STATE_FIPS SQUARE_MIL
1 Ashe County 37009 37 429.350
2 Alleghany County 37005 37 236.459
3 Surry County 37171 37 538.863
4 Gates County 37073 37 342.340
5 Currituck County 37053 37 263.871
6 Stokes County 37169 37 455.793
7 Camden County 37029 37 240.615
8 Warren County 37185 37 443.659
9 Northampton County 37131 37 550.581
10 Hertford County 37091 37 355.525
geometry
1 MULTIPOLYGON (((-81.65649 3...
2 MULTIPOLYGON (((-81.30999 3...
3 MULTIPOLYGON (((-80.71416 3...
4 MULTIPOLYGON (((-76.91183 3...
5 MULTIPOLYGON (((-75.82778 3...
6 MULTIPOLYGON (((-80.43315 3...
7 MULTIPOLYGON (((-76.54193 3...
8 MULTIPOLYGON (((-77.91907 3...
9 MULTIPOLYGON (((-77.16403 3...
10 MULTIPOLYGON (((-77.15428 3...
Simple feature collection with 100 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32186 ymin: 33.84175 xmax: -75.46003 ymax: 36.58815
Geodetic CRS: NAD83
# A tibble: 100 × 9
AREA PERIMETER COUNTYP010 STATE COUNTY FIPS
<dbl> <dbl> <dbl> <chr> <chr> <chr>
1 0.112 1.61 1994 NC Ashe C… 37009
2 0.0616 1.35 1996 NC Allegh… 37005
3 0.140 1.77 1998 NC Surry … 37171
4 0.0891 1.43 1999 NC Gates … 37073
5 0.0687 4.43 2000 NC Currit… 37053
6 0.119 1.40 2001 NC Stokes… 37169
7 0.0626 2.11 2002 NC Camden… 37029
8 0.115 1.46 2003 NC Warren… 37185
9 0.143 2.40 2004 NC Northa… 37131
10 0.0925 1.81 2005 NC Hertfo… 37091
# … with 90 more rows, and 3 more variables:
# STATE_FIPS <chr>, SQUARE_MIL <dbl>,
# geometry <MULTIPOLYGON [°]>
sf classessf [100 × 9] (S3: sf/tbl_df/tbl/data.frame)
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA
..- attr(*, "names")= chr [1:8] "AREA" "PERIMETER" "COUNTYP010" "STATE" ...
[1] "sf" "tbl_df" "tbl"
[4] "data.frame"
[1] "sfc_MULTIPOLYGON" "sfc"
[1] "XY" "MULTIPOLYGON" "sfg"
Coordinate Reference System:
User input: NAD83
wkt:
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]]
Coordinate Reference System:
User input: NAD83 / UTM zone 15N
wkt:
PROJCRS["NAD83 / UTM zone 15N",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["UTM zone 15N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-93,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",26915]]
List of 11
$ : int(0)
$ : int(0)
$ : int(0)
$ : int(0)
$ : int(0)
$ : int 268
$ : int 717
$ : int(0)
$ : int(0)
$ : int(0)
$ : int(0)
- attr(*, "predicate")= chr "intersects"
- attr(*, "region.id")= chr [1:11] "1" "2" "3" "4" ...
- attr(*, "remove_self")= logi FALSE
- attr(*, "retain_unique")= logi FALSE
- attr(*, "ncol")= int 940
- attr(*, "class")= chr [1:2] "sgbp" "list"
logi [1:100, 1:940] FALSE FALSE FALSE FALSE FALSE FALSE ...
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[4,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[7,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[8,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[9,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[10,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[11,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[,8] [,9] [,10] [,11]
[1,] FALSE FALSE FALSE FALSE
[2,] FALSE FALSE FALSE FALSE
[3,] FALSE FALSE FALSE FALSE
[4,] FALSE FALSE FALSE FALSE
[5,] FALSE FALSE FALSE FALSE
[6,] FALSE TRUE FALSE FALSE
[7,] FALSE FALSE FALSE FALSE
[8,] FALSE FALSE FALSE FALSE
[9,] FALSE FALSE FALSE FALSE
[10,] FALSE FALSE FALSE FALSE
[11,] FALSE FALSE FALSE FALSE
nc_air = nc %>%
mutate(
n_air = map_int(
st_intersects(nc, air),
length
)
) %>%
filter(n_air > 0)
nc_air %>% pull(COUNTY) [1] "Forsyth County" "Guilford County"
[3] "Dare County" "Wake County"
[5] "Pitt County" "Catawba County"
[7] "Buncombe County" "Wayne County"
[9] "Mecklenburg County" "Moore County"
[11] "Cabarrus County" "Lenoir County"
[13] "Craven County" "Cumberland County"
[15] "Onslow County" "New Hanover County"
air_nc = air %>%
slice(
st_intersects(nc, air) %>%
unlist() %>%
unique()
)
air_nc %>% pull(AIRPT_NAME) [1] "SMITH REYNOLDS AIRPORT"
[2] "PIEDMONT TRIAD INTERNATIONAL AIRPORT"
[3] "DARE COUNTY REGIONAL AIRPORT"
[4] "RALEIGH-DURHAM INTERNATIONAL AIRPORT"
[5] "PITT-GREENVILLE AIRPORT"
[6] "HICKORY REGIONAL AIRPORT"
[7] "ASHEVILLE REGIONAL AIRPORT"
[8] "SEYMOUR JOHNSON AIR FORCE BASE"
[9] "CHARLOTTE/DOUGLAS INTERNATIONAL AIRPORT"
[10] "MOORE COUNTY AIRPORT"
[11] "CONCORD REGIONAL AIRPORT"
[12] "KINSTON REGIONAL JETPORT AT STALLINGS FIELD"
[13] "CHERRY POINT MARINE CORPS AIR STATION /CUNNINGHAM FIELD/"
[14] "COASTAL CAROLINA REGIONAL AIRPORT"
[15] "POPE AIR FORCE BASE"
[16] "FAYETTEVILLE REGIONAL/GRANNIS FIELD"
[17] "ALBERT J ELLIS AIRPORT"
[18] "WILMINGTON INTERNATIONAL AIRPORT"
Simple feature collection with 13 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32187 ymin: 33.84452 xmax: -75.45998 ymax: 36.58812
Geodetic CRS: WGS 84
# A tibble: 13 × 3
ID DISTRICT geom
<chr> <chr> <MULTIPOLYGON [°]>
1 037108112001 1 (((-77.32845 35.35031, -…
2 037108112002 2 (((-78.89928 35.12619, -…
3 037108112003 3 (((-75.68266 35.23291, -…
4 037108112004 4 (((-78.77926 35.78568, -…
5 037108112005 5 (((-79.8968 36.38075, -7…
6 037108112006 6 (((-80.4201 35.68953, -8…
7 037108112007 7 (((-77.59169 34.40907, -…
8 037108112008 8 (((-78.93373 34.95909, -…
9 037108112009 9 (((-80.93058 35.18181, -…
10 037108112010 10 (((-81.04032 35.40447, -…
11 037108112011 11 (((-84.00768 35.37262, -…
12 037108112012 12 (((-80.4996 35.6493, -80…
13 037108112013 13 (((-79.90775 36.3818, -7…
The Reock score is a measure of compactness that is calculated as the ratio of a shape’s area to the area of its minimum bounding circle.
Simple feature collection with 13 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 196.7724 ymin: 3748.368 xmax: 1002.17 ymax: 4057.317
CRS: +proj=utm +zone=17 +datum=NAD83 +units=km +no_defs
# A tibble: 13 × 4
ID DISTR…¹ geom reock
<chr> <chr> <MULTIPOLYGON [km]> [1]
1 037108… 12 (((545.2987 3945.168, 54… 0.116
2 037108… 13 (((597.9665 4026.852, 59… 0.237
3 037108… 3 (((984.0709 3911.853, 98… 0.266
4 037108… 2 (((691.4141 3889.056, 69… 0.303
5 037108… 9 (((506.3207 3893.208, 50… 0.339
6 037108… 8 (((688.6584 3870.456, 68… 0.342
7 037108… 11 (((226.7522 3918.52, 226… 0.344
8 037108… 6 (((552.4691 3949.669, 55… 0.378
9 037108… 1 (((833.677 3918.083, 831… 0.378
10 037108… 5 (((598.9504 4026.746, 59… 0.399
11 037108… 10 (((496.339 3917.9, 472.7… 0.411
12 037108… 4 (((700.7063 3962.453, 70… 0.480
13 037108… 7 (((813.3004 3812.784, 81… 0.624
# … with abbreviated variable name ¹DISTRICT
Sta 523 - Fall 2022